DRAFT

Ajitsuke Tamago - The art of boiling eggs

Steven Fowler (see more projects here)

No description has been provided for this image

Ajitsuke Tamago, aka Ramen eggs, are one of my favorite additions to a hearty ramen bowl. These eggs have a thick treacly yoke with a soft custard like white that has been marinated. The hardest part of making these eggs is the cooking to get the right egg consistencies.

The goal of this investigation is to model the cooking step of making Ajitsuke Tamago , get an understanding of how the egg undergoes the phase changes that occur. The model uses finite-differences as a method to implement a solution to a PDE ove a 2D space. The ease and simplicity of this implementation can be useful when wanting quick estimates to input into design activities.

What I did

  • I researched the chemistry of eggs to gain an understanding of the mechanisms that drive their phase changed.
    • From this research I defined 3 requirements needed to cook a Ajitsuke Tamago
  • Based on models in the literature, I wrote a function that creates a mesh
    • This mesh defines that geometry of the egg and the varying diffusion coefficients of its constituent parts.
  • The cooking of Ajitsuke Tamago was then simulated:
    • Initial investigation looked into 'over cooking' the egg to get an general idea of the heat transfer response
    • Then simulations with varying rates of heating were undertaken to see if what impact this had on the established requirements

Why I did it

Apart from wanting to explore my love of ramen and Ajitsuke Tamago, I've always been interested in the effort vs gain of modelling and simulation. So, I used this investigation to show how a simple technique can be used to gain meaningful insights.

What I learnt

  • The use of finite-differences to model a diffusion problem gives results that are comparable to the general cooking procedure and other models found in the literature.
  • The variation in heating rates (boiling water around the egg) only impacts the cooking time.
    • It's assumed (not modelled) that very slow rates of heating or effectively lower maximum temperatures would impact the cooking profile. This is seen with sous vide eggs11.
  • The transition of the egg from boiling water to ice water when the yolk-white interface reaches 63°C create a egg that meets the requirements.
    • This takes 7 minutes of boiling times for egg size investigated, which matches the general cooking procedure for Ajitsuke Tamago.

The chemistry of eggs

There are many publications on eggs, their chemistry and even modelling the boiling of eggs 1, 2. In summary, the composition and transformation of cooking an egg is as follows:

Based on Hen eggs, the main areas of an egg consists of 3 key parts:

  1. Shell - This is a porous structure mainly consisting of calcium carbonate and it accounts of 12% of an eggs mass.
  2. White (albumen) - Multi-layered opalescent substance mainly consisting of 2 proteins (12% ovotransferrin, 54% ovalbumin). It accounts for 67% of the egg's mass.
  3. Yoke - Contains fats, vitamins and < 50% of the egg's protein. It accounts for 33% of the egg's mass.

Other parts:

  • Inner and outer membrane between shell and albumen.
  • Chalazae - anchors the yoke in the centre of the egg.
  • Vitelline membrane - between the yolk and albumen.
  • Air cell - Generally forms at the larger end of the egg between the inner and outer membranes. It forms when the egg contracts after being laid.

How the parts of an egg cook (or not)

The phase change of an egg (yolk and white) is due to a coagulation of denatured proteins within the egg. The proteins denature at different temperature ranges, aggregate and form gel like structures. Note: Proteins can also be denatured by other mechanisms (i.e. pH changes, mechanical action) but these are out-of-scope of this investigation.

The phase changes of an egg are summarized in the below table (based on ref 3):

Temperature (°C) White Yolk
63 Partially set Runny thickening
65 - 70 Soft set gel Runny treacle
73 Largely set Soft gel
77 Tender solid Hard boiled
80 Firm Onset green discoloration
90 Over cooked Dry, crumbly solid

A general recipe for Ajitsuke Tamago is as follows4:

  1. Boil water (100°C)
  2. Place egg in water and cook for 7 minutes
  3. Immerse egg in ice water (0°C)
  4. Once cooled remove shells
  5. Marinade eggs for 24 hrs

Modelling requirements

Based on the above understanding, an Ajitsuke Tamago requires:

ID_req Requirement Criteria
1 Yolk temperature shall not exceed 70°C to keep yolk runny max(yolk temperature) <= 70°C
2 Yolk temperature shall not exceed 65°C to thicken the yolk max(yolk temperature) <= 63°C
3 White temperature shall not overly (on average) exceed 80°C to prevent onset of discoloration mean(white temperature) =< 80°C

Modelling Ajitsuke Tamago

The section provides an overview of the model used to simulate the Ajitsuke Tamago cooking. There are more detailed notes in the appendix. In summary, the following is covered:

  • An overview of the maths deriving the model
  • The creation of the mesh for the model
    • Defining the egg shape and areas of interest
    • Creating the mesh (diffusion coefficient array)
  • The implementation of the model using the mesh

Overview of modelling math

To model the cooking of an egg, a general diffusion model was implemented to simulate the heat transfer in the system. A 2D diffusion type system has the general form of: $$ \frac{\delta u}{\delta t} = \alpha (\frac{\delta^2 u}{\delta x^2} + \frac{\delta^2 u}{\delta i^2}) $$ where $u(x,i,t)$ is the function to be solved, $j$ and $i$ are cartesian coordinates and t is time. $\alpha$ is the diffustion coefficient and determines the rate change in $u$ in $t$.

Based on finite differences (see appendix for detailed notes), this can be expressed as: $$ u_i^{n+1} = u_{i,j}^n + F_x(u_{i+1,j}^n - 2u_{i,j}^n + u_{i-1,j}^n) + F_y(u_{i,j+1}^n - 2u_{i,j}^n + u_{i,j-1}^n) $$ where $$ F_x = \alpha \frac{\Delta t}{\Delta j^2}, F_y = \alpha \frac{\Delta t}{\Delta i^2} $$

$F$ is the Fourier Number, a dimensionless value that represents the physical parameters of the diffusion problem in $\alpha$ and the discretized parameters $\Delta j$ and $\Delta t$.

The time domain is defined as: $$ t_n = n \Delta t, n = 0,...,N_t $$ with the space domain as: $$ x_i = i\Delta x, i = 0,...,N_x, $$ $$ y_j = j\Delta i, j = 0,...,N_y $$

Modelling Mesh

The standard mesh for a finite difference problem is a cartesian grid. Although it is possible to make more complex or dynamic meshes, assuming this is needed to capture more detail in an simulation then an FEA approach is probably more desirable. Also, it is acknowledged that the symmetry of the model means mesh size and computations could be reduced, however this will be seen as a future optimisiation. Assuming that the relatively high diffusion coefficient (≈ x10 greater than yolk and white) and a thickness of ≈0.3 mm of the shell has little its influence on the overall heat transfer, it can be ignored. So, the overall model consists of 2 parts: the white and the yolk.

Egg shape

The shape of the egg based on the model from Lian et al5, where the parametric form is: $$ y_p=asin\theta $$ $$ x_p=bcos\theta(1+c_1sin\theta+c_2sin^2\theta+c_3sin^3\theta) $$ where $y_p$ and $x_p$ are the coordinates of the egg profile (2D); $\theta$ is the eccentric angle; $a$ the half length; $b$ approximate half width; $c_1,c_2,c_3$ are parameters to be adjusted to fit an egg shape.

Show code
In [1]:
# Egg shape curve
# A parametric equation to generate an egg shape curve
# Reference: Comparison of egg-shape equations using relative curvature measures of nonlinearity, Lian et al., 2024

import numpy as np
import matplotlib.pyplot as plt

# Parametic equations that define egg shape
#   - Equation of ellipse with 3rd order parameters
def egg_shape(half_width, half_height, theta, c1, c2, c3 , x_axis_offset=0, y_axis_offset=0):
    """Parametric equation that defines the egg shape
    Function is that of an ellipse with 3rd order parameters that can be used to fit egg shape

    Args:
        half_width (float): width of egg, radii
        half_height (float): height of egg, radii
        theta (float): eccentric angle, theta
        c1 (_type_): 1st order parameter of egg fitting
        c2 (_type_): 2nd order parameter of egg fitting
        c3 (_type_): 3rd order parameter of egg fitting
        x_axis_offset (int, optional): Shift in x-axis. Defaults to 0.
        y_axis_offset (int, optional): Shift in y-axis. Defaults to 0.

    Returns:
        tuple: x, y - coordinates of point on egg shape corresponding to theta.
    """
    x = (half_width * np.cos(theta)) * (1 + c1*np.sin(theta) + c2*np.sin(theta)**2 + c3*np.sin(theta)**3) + x_axis_offset
    y = half_height * np.sin(theta) + y_axis_offset
    return x, y

# Returns the difference (error) between the magnitude of a point (point_x, point_y) and the corresponding 
# point on the egg shape.
# The corresponding point is that which has the same angle theta. 
def find_error_at_t(point_x, point_y, half_width, half_height, c1 , c2 , c3 , t, x_axis_offset, y_axis_offset):
    egg_x, egg_y = egg_shape(half_width, half_height, t, c1, c2, c3, 0, 0)
    err = np.atan2(point_y - y_axis_offset, point_x - x_axis_offset) - np.atan2(egg_y, egg_x)
    return err, egg_x, egg_y

# Use of bisection method to find the corresponding point on the egg shape. Where the corresponding point 
# is that which has the same angle theta. 
def find_theta_of_egg(point_x, point_y, half_width, half_height, c1 , c2 , c3 , x_axis_offset, y_axis_offset):
    a = 0
    b = np.pi 
    err_a = find_error_at_t(point_x, point_y, half_width, half_height, c1 , c2 , c3 , a, x_axis_offset, y_axis_offset)[0]
    err_b = find_error_at_t(point_x, point_y, half_width, half_height, c1 , c2 , c3 , b, x_axis_offset, y_axis_offset)[0]
    
    # Ensuring guess a & b are either side of the solution
    if np.sign(err_a)  == np.sign(err_b):
        b = -np.pi
        err_b = find_error_at_t(point_x, point_y, half_width, half_height, c1 , c2 , c3 , b, x_axis_offset, y_axis_offset)[0]
        if np.sign(err_a)  == np.sign(err_b):
                raise Exception("Bisection method requires a sign change over the interval.")
    
    # Bisection method
    count = 0
    while True:
        m = (a + b) / 2
        err_m, egg_x, egg_y = find_error_at_t(point_x, point_y, half_width, half_height, c1 , c2 , c3 , m, x_axis_offset, y_axis_offset)

        if np.abs(err_m) < 1e-14:
            return m, egg_x, egg_y
        elif np.sign(err_a) == np.sign(err_m):
            a = m
            err_a = err_m
        else:
            b = m
            err_b = err_m
        
        count += 1
        if count >= 1000:
            print("Error: ", err_m)
            raise Exception("Bisection method did not converge within 1000 iterations.")

# Find if a point (x, y) is equal to or within the parameter of the egg shape
# Return boolean
def point_within_egg(point_x, point_y, half_width, half_height, c1 , c2 , c3 , x_axis_offset, y_axis_offset):

    # Find a theta of the egg shape that matches the theta of the point
    egg_theta, egg_x, egg_y = find_theta_of_egg(point_x=point_x, point_y=point_y, half_width=half_width, half_height=half_height,
                                                c1=c1, c2=c2, c3=c3, x_axis_offset=x_axis_offset, y_axis_offset=y_axis_offset)
    
    # Return if point is within egg shape
    return ((point_x - x_axis_offset)**2 + (point_y - y_axis_offset)**2) <= (egg_x**2 + egg_y**2), egg_x, egg_y

# Find if a point (x, y) is equal to or within the parameter of the circle
# Assumes circle is the shape of the yolk
# Return boolean
def yoke_shape(x_axis_offset, y_axis_offset, a, b, t):
    x = x_axis_offset + a*np.cos(t)
    i = y_axis_offset + b*np.sin(t)
    return x, i

#######################################################################
#                                                                     #
# Testing egg shape functions and allocation of points defined above  #
#                                                                     #
#######################################################################

t = np.linspace(0, 2*np.pi, 100)
# Egg parameters
width = 15
height = 25
c1 = -0.2
c2 = 0.1
c3 = 0.0
x_shift = width
y_shift = height
# Yoke parameters
a = 10
b = 10
x_shift_yoke = x_shift
y_shift_yoke = y_shift - 10

# Egg and yolk perimeters
x, y = egg_shape(width, height, t, c1=c1, c2=c2, c3=c3, x_axis_offset=x_shift, y_axis_offset=y_shift)
x_yoke, y_yoke = yoke_shape(x_shift_yoke, y_shift_yoke, a, b, t)

# Random test points
x_test = np.random.rand(1000) * width * 3 - 0.5 * width
y_test = np.random.rand(1000) * height * 3 - 0.5 * height

# Plotting
point_in_egg = np.zeros_like(x_test, dtype=bool)
point_in_york = np.zeros_like(x_test, dtype=bool)
x_egg_chk = np.zeros_like(x_test)
y_egg_chk = np.zeros_like(x_test)
for i in range(len(x_test)):
    point_in_egg[i], x_egg_chk[i], y_egg_chk[i] = point_within_egg(x_test[i], y_test[i], width, height, 
                                       c1 , c2 , c3 , 
                                       x_shift, y_shift)
    point_in_york[i] = ((x_test[i] - x_shift_yoke)**2 + (y_test[i] - y_shift_yoke)**2)**0.5 <= (a)
    

plt.plot(x_test[np.where(np.invert(point_in_egg))], y_test[np.where(np.invert(point_in_egg))], 'go')
plt.plot(x_test[np.where(point_in_egg != point_in_york)], y_test[np.where(point_in_egg != point_in_york)], 'ro')
plt.plot(x_test[np.where(point_in_york)], y_test[np.where(point_in_york)], 'yo') 
plt.plot(x_egg_chk, y_egg_chk, 'bo')
plt.plot(x, y, 'k-')
plt.plot(x_yoke, y_yoke, 'k-')
plt.title('Egg Shape Curve - simulation using random points')
plt.xlabel('Width (mm)')
plt.ylabel('Height (mm)')
plt.axis('equal')
plt.style.use('seaborn-v0_8-whitegrid')
plt.show()
No description has been provided for this image

Figure 1. Illustration of egg shape fitting function.
The blue line shows the egg like shape that has been generated by the function. Randomly generated points are allocated to the specific area of the array/image: Yellow points are within the yolk; Red points are within the white and; green points are outside the egg shape. The blue points are those on the same vector as the randomly generated points (plotted from the origin).

The above plot shows the egg shape using fictional parameters to simulate an egg. The points in the plot are randomly generated and color-assigned depending on where are in the egg. This approach will be used to generate the mesh and thermal diffusion coefficient array for the modelling.

Define thermal diffusion coefficient array

The model will be defined by an array of varying diffusion coefficients that represent each area of the egg (white and yolk). Using the functions defined in the above section, the parameters of the egg fitting function were adjusted to fit measured data (see table below) and the yolk size was estimated based on the mass of the white and yolk (yolk being ≈33% of white + yolk).

Three measures were sampled on medium sized supermarket eggs. The parameters of the egg shape equation were then adjusted to fit these measurements.

Measurement Value (mm)
Max width 45.8
Heigh of 'max width' 25
Max height 55.8

The thermal properties of the egg were found in publications (see table below). However, as can be seen in the table the quoted magnitude of these properties varies considerably.

Egg part Reference Density (Kg.cm-3) Conductivity (W.m-1.K-1) Capacity (J.Kg-1.K-1)
Yolk 6 1028 0.54 3700
Yolk 7 1100 0.31 3700
Yolk* 8 1036 0.34 2093
Yolk* 8 1036 0.34 3266
Yolk 9 N/A 0.34 3180
Yolk 10 1032 0.34 2700
White 7 1134 0.48 3800
White* 8 1036 0.59 2093
White* 8 1036 0.59 3935
White 9 N/A 0.58 3180
White 10 1038 0.54 3700
  • Reference [8] quotes a range of values, hence there are 2 rows for each egg part.

For the purposes of this model, the values from reference [11] were chosen:

Egg part Density (Kg.cm-3) Conductivity (W.m-1.K-1) Capacity (J.Kg-1.K-1) Diffusion (mm2.s-1)
Yolk 1032 0.34 2700 0.12
White 1038 0.54 3700 0.14
Show code
In [2]:
import math
import numpy as np

# Rounds value to nearest odd number
# Returns int.
# Not used in method
def round_up_to_nearest_odd(value):
    round_up_value = math.ceil(value)
    if round_up_value % 2 == 0:
        return int(round_up_value +1)
    else:
        return int(round_up_value)
    
# Model definition
egg_width = 45.8                    # mm
egg_height = 55.8                   # mm
egg_c1 = -0.1                       # Shape parameter 1st order
egg_c2 = 0.1                        # Shape parameter 2nd order
egg_c3 = 0.0                        # Shape parameter 3rd order
white_alpha = 0.14060303077644118   # Egg white coefficient of diffusion
yoke_width = egg_width * 0.41**0.5  # mm, estimate of 33% by mass (area)
yoke_height = yoke_width            # mm, assume round for now 
yoke_center = -5                    # mm, height from base off egg
yoke_alpha  =0.1220212460522538     # Egg yoke coefficient of diffusion

# Mesh definition
dx = 0.5     # mm per node
dy = dx

mesh_X_mm = np.linspace(-egg_width/2 - 1, egg_width/2 + 1, int(np.ceil((egg_width + 2)/dx)))
mesh_Y_mm = np.linspace(-egg_height/2 - 1, egg_height/2 + 1, int(np.ceil((egg_height + 2)/dx)))

# Diffusion grid
X, Y = np.meshgrid(mesh_X_mm, mesh_Y_mm, indexing='ij')
alpha_grid = np.zeros_like(X)
boundary_coord = []
white_coord = []
yolk_coord = []

# Allocate each node a diffusion coefficient
for i in range(np.size(alpha_grid,0 )):
    for j in range(np.size(alpha_grid, 1)):
        # Check if point in yoke
        point_in_yoke = (X[i, j]**2 + (Y[i, j] - yoke_center)**2) <= (yoke_width/2)**2
        if point_in_yoke:
            yolk_coord.append((i, j))
        else:
            # Check if point in white
            point_in_egg, e1, e2 = point_within_egg(X[i, j], Y[i, j], egg_width/2, egg_height/2,
                                                    egg_c1, egg_c2, egg_c3, 
                                                    0, 0)
            if point_in_egg:
                white_coord.append((i, j))
            else:
                # Point outside of egg shape
                boundary_coord.append((i, j))

boundary_coord = np.array(boundary_coord)
white_coord = np.array(white_coord)
yolk_coord = np.array(yolk_coord)

print(f"yolk to white ratio: {np.size(yolk_coord,0) / (np.size(white_coord,0) + np.size(yolk_coord,0)):.2f}")

#######################################################################
#                                                                     #
# Plotting mesh                                                       #
#                                                                     #
#######################################################################

i, j = egg_shape(egg_width/2, egg_height/2,
                 np.linspace(0, 2*np.pi, 100),
                 egg_c1, egg_c2, egg_c3,
                 0, 0)

alpha_grid[yolk_coord[:, 0], yolk_coord[:, 1]] = yoke_alpha
alpha_grid[white_coord[:, 0], white_coord[:, 1]] = white_alpha
alpha_grid[boundary_coord[:, 0], boundary_coord[:, 1]] = 0

fig, ax = plt.subplots()
ax.set_aspect('equal')
scatter = ax.scatter(X, Y, c=alpha_grid, cmap='viridis', s=0.1)
fig.colorbar(scatter, ax=ax, label='Diffusion coefficient')
ax.plot(i, j)
ax.set_title('Image of mesh (thermal diffusion coefficient array)')
ax.set_xlabel('X (mm)')
ax.set_ylabel('Y (mm)')
plt.style.use('seaborn-v0_8-whitegrid')
plt.grid(False)
plt.show()

yolk to white ratio: 0.33
No description has been provided for this image

Figure 2. Illustration of diffusion coefficient array.
This plot shows the nodes of the cartesian mesh (scaled to mm), where the colour of the node represents the magnitude of the diffusion coefficient. The blue line highlights the boundary of the egg. Note: Nodes with 0 diffusion coefficient (dark purple) are the boundary nodes.

The above image illustrates the diffusion array for the model. Outside of the egg is the boundary nodes with a coefficient of 0, the yellow (≈0.14 mm2s-1) is the white and green (≈0.12 mm2s-1) is the yolk. The origin is at the centre of the egg-shape function and there is a 0.5 mm resolution (node spacing).

Finite-difference model

A generic model for finite differences was drafted (see appendix) for solving general diffusion problem. However, for this investigation an diffusion coefficient array was needed to define the circular egg shape (without using a polar coordinate system) and dynamic control was needed to investigate the requirements. Hence, the below ajitsuke_egg_simulation function was written.

Criteria needed to investigate requirements:

  • Geometrically defined diffusion coefficients
  • Dynamic heating (control of boundary conditions)
  • Max temp monitoring (logic to drive heating temperature)
  • Mean yolk monitoring (logic to drive heating temperature)

Model assumptions

ID_ma Assumption
1 Diffusion coefficients are constant throughout process (see appendix - this is not true)
2 The error in the finite-difference approach is within the noise of the system
3 The egg shape defined in the coefficient array is representative of a general egg (small samples numbers measured)
4 A uniform boundary condition will give representative results i.e. ignoring the dynamics of the boiling water
5 Shell and other egg parts (see the chemistry of eggs section above) have a minimal impact of the thermal diffusion within the egg
Show code
In [3]:
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.cm as cm
import time as time
from IPython.display import HTML
import numpy as np

def ajitsuke_egg_simulation(dt, nt, dx, dy, alpha_grid, 
                            heat_rate, mean_yolk_temp_switch, max_yolk_temp_switch,
                            boundary_T_cool, boundary_T_max, boundary_T_initial,
                            egg_T_initial, boundary_coord, yolk_coord, white_coord):
    """Use of finite-differences to solve diffusion problem with 2D coefficient array input

    Args:
        dt (float): Step size in time. 
        nt (int): Number of steps in time.
        dx (float): Step size in x-axis
        dy (float): Step size in y-axis
        alpha_grid (arr): Array of diffusion coefficients. Array size defines mesh.
        heat_rate (float): Rate of boundary temperature increase from t0 to boundary_T_max. 
        mean_yolk_temp_switch (float): Mean yolk temperature that boundary conditions change at.
        max_yolk_temp_switch (float): Max yolk (yolk-white interface) temperature that boundary conditions change at.
        boundary_T_cool (float): Temperature of boundary conditions for cooling.
        boundary_T_max (float): Maximum temperature of boundary conditions.
        boundary_T_initial (float): Initial temperature of boundary conditions.
        egg_T_initial (float): Initial egg temperature. Assumes uniform initial temperature over all of egg.
        boundary_coord (lst): Coordinate list of boundary nodes.
        yolk_coord (lst): Coordinate list of yolk nodes.
        white_coord (lst): Coordinate list of white nodes.

    Raises:
        ValueError: If stability conditions are not met.

    Returns:
        arr: 3D array. u(t, i, j).
    """

    start_time = time.time()

    nx = np.size(alpha_grid, 0)
    ny = np.size(alpha_grid, 1)

    # Check stability condition - using max alpha for worst case
    if dt > (dx**2 * dy**2) / (2 * alpha_grid.max() * (dx**2 + dy**2)):
        dt_max = (dx**2 * dy**2) / (2 * alpha_grid.max() * (dx**2 + dy**2))
        alpha_max = (dx**2 * dy**2) / (2 * dt * (dx**2 + dy**2))
        raise ValueError("Time step size dt is too large for stability. Maximum dt: ", dt_max, " or reduce alpha to: ", alpha_max)


    # Set initial conditions
    boundary_T_current = boundary_T_initial
    u0 = np.empty_like(alpha_grid) 
    u0[boundary_coord[:, 0], boundary_coord[:, 1]] = boundary_T_current
    u0[white_coord[:, 0], white_coord[:, 1]] = egg_T_initial
    u0[yolk_coord[:, 0], yolk_coord[:, 1]] = egg_T_initial

    # Boundary conditions with time
    cool = True
    mean_yolk = None
    mean_white = None
    u = np.empty((nt, nx, ny))  # Initialize the solution array
    u[0, :, :] = u0  # Set initial condition
    for n in range(1, nt):
        
        mean_yolk = np.mean(u[n-1, yolk_coord[:, 0], yolk_coord[:, 1]])
        max_yolk = np.max(u[n-1, yolk_coord[:, 0], yolk_coord[:, 1]])
        mean_white = np.mean(u[n-1, white_coord[:, 0], white_coord[:, 1]])
        # Real time data print
        # print(f"Mean yolk: {mean_yolk:.2f}, Mean white: {mean_white:.2f} - after {dt*n:.2f} seconds\r", end='')
        
        # Cut off heat source when requirement 2 is met (mean_yolk_temp_switch)
        if ((mean_yolk >= mean_yolk_temp_switch) or (max_yolk >= max_yolk_temp_switch)) and (cool):
            boundary_T_current = boundary_T_cool
            u[n-1, boundary_coord[:, 0], boundary_coord[:, 1]] = boundary_T_current
            cool = False
            print(f"Temperature changed! - Mean yolk: {mean_yolk:.2f}, Mean white: {mean_white:.2f} - after {dt*n:.2f} seconds")

        # Cap heating to boundary_T_max
        if (boundary_T_current < boundary_T_max) & (cool):
            boundary_T_current += heat_rate
            u[n-1, boundary_coord[:, 0], boundary_coord[:, 1]] = boundary_T_current

        # Boundary condition check
        u[n, boundary_coord[:, 0], boundary_coord[:, 1]] = boundary_T_current

        # Iterate over internal grid points
        for i in range(1, nx-1):    
            for j in range(1, ny-1):
                u[n, i, j] = (u[n-1, i, j] + 
                            alpha_grid[i, j] * dt * ((u[n-1, i+1, j] - 2*u[n-1, i, j] + u[n-1, i-1, j]) / dx**2 +
                                                    (u[n-1, i, j+1] - 2*u[n-1, i, j] + u[n-1, i, j-1]) / dy**2))
                
    end_time = time.time()
    elapsed_time = end_time - start_time
    # print(f"Mean yolk: {mean_yolk:.2f}, Mean white: {mean_white:.2f} - after {dt*n:.2f} seconds")
    print(f"Elapsed time for solution: {elapsed_time:.2f} seconds.\n Time simulated: {dt * nt:.2f} seconds.")

    return u

Simulation

The following simulation initially looks at constant heating of the egg to an over cooked state. Then different heating regimes are simulated with the egg being cooled in an ice bath when the yolk-white interface reached 63°C.

Over cooking - constant heating

As an initial look at the model, heating of the egg from room temperature (20°C) with an constant boundary condition of 100°C (boiling water) will be simulated.

Show code
In [16]:
dt = 0.035
nt = 35000 # >20 minutes

mean_yolk_temp_switch = 1000    # Too high to trigger boundary condition change
max_yolk_temp_switch = 1000     # Too high to trigger boundary condition change
boundary_T_cool = 0
boundary_T_max = 100
egg_T_initial = 20
heat_rate = 0 * dt
boundary_T_initial = 100

data_100 = ajitsuke_egg_simulation(dt, nt, dx, dy, alpha_grid, 
                                heat_rate, mean_yolk_temp_switch, max_yolk_temp_switch, 
                                boundary_T_cool, boundary_T_max, boundary_T_initial,
                                egg_T_initial, boundary_coord, yolk_coord, white_coord)

# Save for later access without needing to re-run simulation
np.save("egg_100_heating", data_100)

Elapsed time for solution: 599.03 seconds.
 Time simulated: 1225.00 seconds.
Show code
In [4]:
from matplotlib import collections  as mc

# Option to load data rather than re-running simulation
# dt = 0.035
# data_100 = np.load("egg_100_heating.npy")

# Plot data on axes from subplots
# Return x & y data
def plot_simulation_data(data, axes, boundary_coord, white_coord, yolk_coord):
    x = np.zeros(np.size(data, 0))
    y_label = ["boundary T", "mean white T", "mean yolk T", "min yolk T", "max yolk T", "max_white_T"]
    y_colour = ["g-", "y-", "b-", "b--", "b--", "w "]
    y = np.zeros((np.size(y_label), np.size(x)))

    for n in range(np.size(x)):
        x[n] = n * dt
        y[0, n] = np.mean(data[n, boundary_coord[:, 0], boundary_coord[:, 1]])
        y[1, n] = np.mean(data[n, white_coord[:, 0], white_coord[:, 1]])
        y[2, n] = np.mean(data[n, yolk_coord[:, 0], yolk_coord[:, 1]])
        y[3, n] = np.min(data[n, yolk_coord[:, 0], yolk_coord[:, 1]])
        y[4, n] = np.max(data[n, yolk_coord[:, 0], yolk_coord[:, 1]])
        y[5, n] = np.max(data[n, white_coord[:, 0], white_coord[:, 1]])

    for n in range(np.size(y_label)):
        axes.plot(x, y[n], y_colour[n], label=y_label[n])

    axes.fill_between(x, y[3], y[4], facecolor='grey', alpha=0.5)

    axes.hlines([63, 80], xmin=np.min(x), xmax=np.max(x), colors='k', linestyles='-.')
    axes.text(0, 63, '63°C', verticalalignment='bottom', horizontalalignment='left', color='k')
    axes.text(0, 80, '80°C', verticalalignment='bottom', horizontalalignment='left', color='k')
    
    axes.set_title('Water at 100°C from t$_0$')
    axes.set_xlabel('Time (s)')
    axes.set_ylabel('Temperature (°C)')

    axes.grid(False)

    return x, y

plt.style.use('seaborn-v0_8-whitegrid')
fig, ax = plt.subplots()

x, y_100 = plot_simulation_data(data_100, ax, boundary_coord, white_coord, yolk_coord)

fig.suptitle('Ahitsuke Tamago simulation:')
plt.show()
No description has been provided for this image

Figure 3. Simulation of boiling egg at 100°C for > 1200 seconds.
Green is boundary condition (temperature of water cooking the egg); Yellow is the mean temperature of the white; Blue is the mean (solid), max and min (dashed) yolk temperature. Grey area highlights temperature range of yolk. Black dash-dot lines are at: 63°C and 80°C.

Show code
In [5]:
# Useful numbers for discussion below...
time_reach_63 = x[np.argmax(y_100[4] >= 63)]

print("Useful numbers...\n")
print(f"\tx Yolk starts to heat (T > 21°C) at: {x[np.argmax(y_100[4] > 21)]:.2f} seconds")
print(f"\tx Yolk-white interface reaches 63°C: {time_reach_63:.2f} seconds ({time_reach_63/60:.1f} minutes)")
print(f"\tx Mean Yolk temperature when yolk-white interface at 63°C: {y_100[2, np.argmax(x >= time_reach_63)]:.2f}°C")
print(f"\tx Mean White temperature when yolk-white interface at 63°C: {y_100[1, np.argmax(x >= time_reach_63)]:.2f}°C")
print(f"\tx Max White temperature when yolk-white interface at 63°C: {y_100[5, np.argmax(x >= time_reach_63)]:.2f}°C")

Useful numbers...

	x Yolk starts to heat (T > 21°C) at: 38.64 seconds
	x Yolk-white interface reaches 63°C: 421.30 seconds (7.0 minutes)
	x Mean Yolk temperature when yolk-white interface at 63°C: 41.97°C
	x Mean White temperature when yolk-white interface at 63°C: 76.90°C
	x Max White temperature when yolk-white interface at 63°C: 99.22°C

From the above plot it can be seen that the yolk does not start to heat for 39 seconds, it them reaches 63°C in 7 minutes. This is a nice alignment as the recipe for cooking the Ajitsuke Tamago stated a 7 minutes cooking time and reference [6] & [10] defines it's egg simulation complete when the yolk-white interface reaches 63°C. When the yolk-white interface reaches 63°C the mean white temperature is at 77°C, which is below at the 'tender' hardness (firm is at 80°C).

Varied heating rate and cooling at 63°C inflection point

As discussed above, reference [6] [10] assumed that the egg is cooked when the yolk-white interface reached 63°C. Consequently, for the next simulations when this condition is met the boundary condition will be dropped to 0°C simulated the egg being submerged in ice water.

3 differing heating scenarios will be modelled:

  • Scenario 1 - Starting with boiling water (reference)
  • Scenario 2 - Starting at 20°C and heating to 100°C at a rate of 0.50°C.s-1
  • Scenario 3 - Starting at 20°C and heating to 100°C at a rate of 0.25°C.s-1

To simulate the egg being present as the water is being boiled.

Show code
In [20]:
dt = 0.035
nt = 35000 # >20 minutes
 # Defines dimensions of problem
heat_rate = 0.5 * dt
mean_yolk_temp_switch = 100  # Requirement 2
max_yolk_temp_switch = 63
boundary_T_cool = 0
boundary_T_max = 100
boundary_T_initial = 20
egg_T_initial = 20

###
data_050 = ajitsuke_egg_simulation(dt, nt, dx, dy, alpha_grid, 
                                heat_rate, mean_yolk_temp_switch, max_yolk_temp_switch, 
                                boundary_T_cool, boundary_T_max, boundary_T_initial,
                                egg_T_initial, boundary_coord, yolk_coord, white_coord)
# Save for later access without needing to re-run simulation
np.save("egg_050_heating", data_050)

###
heat_rate = 0.25 * dt
data_025 = ajitsuke_egg_simulation(dt, nt, dx, dy, alpha_grid, 
                                heat_rate, mean_yolk_temp_switch, max_yolk_temp_switch,
                                boundary_T_cool, boundary_T_max, boundary_T_initial,
                                egg_T_initial, boundary_coord, yolk_coord, white_coord)
# Save for later access without needing to re-run simulation
np.save("egg_025_heating", data_025)

###
heat_rate = 0 * dt
boundary_T_initial = 100
data_000 = ajitsuke_egg_simulation(dt, nt, dx, dy, alpha_grid, 
                                heat_rate, mean_yolk_temp_switch, max_yolk_temp_switch, 
                                boundary_T_cool, boundary_T_max, boundary_T_initial,
                                egg_T_initial, boundary_coord, yolk_coord, white_coord)
# Save for later access without needing to re-run simulation
np.save("egg_000_heating", data_000)

Temperature changed! - Mean yolk: 42.09, Mean white: 76.91 - after 504.00 seconds
Elapsed time for solution: 605.89 seconds.
 Time simulated: 1225.00 seconds.
Temperature changed! - Mean yolk: 42.46, Mean white: 76.89 - after 592.38 seconds
Elapsed time for solution: 607.22 seconds.
 Time simulated: 1225.00 seconds.
Temperature changed! - Mean yolk: 41.97, Mean white: 76.90 - after 421.33 seconds
Elapsed time for solution: 605.79 seconds.
 Time simulated: 1225.00 seconds.
Show code
In [6]:
# Option to load data rather than re-running simulation
# dt = 0.035
# data_000 = np.load("egg_000_heating.npy")
# data_050 = np.load("egg_050_heating.npy")
# data_025 = np.load("egg_025_heating.npy")
 
plt.style.use('seaborn-v0_8-whitegrid')
fig, axs = plt.subplots(1, 3, figsize=(12, 5), sharey=True)

x, y_000 = plot_simulation_data(data_000, axs[0], boundary_coord, white_coord, yolk_coord)
x, y_050 = plot_simulation_data(data_050, axs[1], boundary_coord, white_coord, yolk_coord)
x, y_025 = plot_simulation_data(data_025, axs[2], boundary_coord, white_coord, yolk_coord)

fig.suptitle('Ahitsuke Tamago simulation:')

plt.show()
No description has been provided for this image

Figure 4. Ahitsuke Tamago cooking simulation with differing boundary conditions.
Green is boundary condition (temperature of water cooking the egg); Yellow is the mean temperature of the white; Blue is the mean (solid), max and min (dashed) yolk temperature. Grey area highlights temperature range of yolk. Black dash-dot lines are at 65°C and 70°C, the limits for the yolk.

Show code
In [7]:
# Useful numbers for discussion below...

print("useful numbers...\n")
print(f"Maximum yolk temperature for:               \n\tScenario 1 -> {np.max(y_000[4]):.2f}°C; \n\tScenario 2 -> {np.max(y_050[4]):.2f}°C; \n\tScenario 3 -> {np.max(y_025[4]):.2f}°C.\n")
print(f"Maximum white temperature for:              \n\tScenario 1 -> {np.max(y_000[5]):.2f}°C; \n\tScenario 2 -> {np.max(y_050[5]):.2f}°C; \n\tScenario 3 -> {np.max(y_025[5]):.2f}°C.\n")
print(f"Maximum mean white temperature reached for: \n\tScenario 1 -> {np.max(y_000[1]):.2f}°C; \n\tScenario 2 -> {np.max(y_050[1]):.2f}°C; \n\tScenario 3 -> {np.max(y_025[1]):.2f}°C.\n")
print(f"Minimum egg (yolk) temperature reached for: \n\tScenario 1 -> {np.max(y_000[3]):.2f}°C; \n\tScenario 2 -> {np.max(y_050[3]):.2f}°C; \n\tScenario 3 -> {np.max(y_025[3]):.2f}°C.\n")
print(f"Time of 63°C inflection point:              \n\tScenario 1 -> {x[np.argmax(y_000[0] == 0)]/60:.1f} min; \n\tScenario 2 -> {x[np.argmax(y_050[0] == 0)]/60:.1f} min; \n\tScenario 3 -> {x[np.argmax(y_025[0] == 0)]/60:.1f} min.\n")

useful numbers...

Maximum yolk temperature for:               
	Scenario 1 -> 64.24°C; 
	Scenario 2 -> 64.25°C; 
	Scenario 3 -> 64.27°C.

Maximum white temperature for:              
	Scenario 1 -> 99.22°C; 
	Scenario 2 -> 99.23°C; 
	Scenario 3 -> 99.22°C.

Maximum mean white temperature reached for: 
	Scenario 1 -> 76.90°C; 
	Scenario 2 -> 76.91°C; 
	Scenario 3 -> 76.89°C.

Minimum egg (yolk) temperature reached for: 
	Scenario 1 -> 39.34°C; 
	Scenario 2 -> 39.45°C; 
	Scenario 3 -> 39.80°C.

Time of 63°C inflection point:              
	Scenario 1 -> 7.0 min; 
	Scenario 2 -> 8.4 min; 
	Scenario 3 -> 9.9 min.

Looking at the plot and reinforced by the 'useful numbers' above, apart from the initial heating of the egg the response and temperature at and after the 63°C inflection are all the same for the 3 scenarios. Consequently, the only difference the length of time until the 63°C inflection points which are at 7.0, 8.4 & 9.9 minutes for Scenario 1, 2, & 3 respectively.

Requirements

Considering the above plots in context of the requirements:

  1. Yolk temperature shall not exceed 70°C to keep yolk runny
    • This was met for all scenarios as the maximum yolk temperature was ≈ 64°C.
  2. Yolk temperature shall exceed 65°C to thicken the yolk
    • This was not met for all scenarios as the maximum temperature was ≈ 64°C.
  3. White temperature shall not overly exceed 80°C to prevent onset of discoloration
    • This was met for all scenarios as the maximum mean white temperature was ≈ 77°C.

In summary, all scenarios produce the same egg. However, it could be argued that scenario 1 is less prone to error as you don't need to factor the rate of water heating into the cooking.

Conclude

The above report demonstrates that the use of finite-differences can provide an quick and informative model to simulation a thermal diffusion problem. The results of the simulation align with standard cooking procedure (cooking time of 7 minutes in boiling water) and the removal of the egg when the yolk-white interface is at 63°C results in the simulation meeting the proposed requirements.

Further work & Improvements (some ideas)

  • Conduct a parameter analysis to establish the 'sensitivity' of input parameters to the perfect egg
  • Once a parameter analysis is complete, there should be data to conclude on the operating window in which an acceptable egg is produced
  • Add an egg fitting algorithm based on images of eggs
    • This would facilitate the increased sampling of egg distributions
  • With the aforementioned ability to measure egg size, include the egg size in the the parameter analysis
  • Data from cooking eggs (lab work!) would be good to help validate the simulation results.
    • This would be particularly good if it included varying egg sizes.

Heat profile

The following are additional plots and videos from the simulations.

Show code
In [ ]:
# Profile of egg at 63°C inflection point
index_reach_63 = np.argmax(y_000[4] >= 63)
crosssection_point = 50

fig, (ax1, ax2) = plt.subplots(1, 2)

# Left hand plot
ax1.text(5, crosssection_point, 'Cross-section', verticalalignment='bottom', horizontalalignment='left', color='b')
ax1.axhline(crosssection_point, linestyle='--', label='plot line', color='blue')
im1 = ax1.imshow(data_000[index_reach_63,].T, cmap='hot', vmin=0, vmax=100, origin='lower') # Transpose needed as array is [i,j] not [x,y]
ax1.set_xlabel('Size (nodes)')
ax1.set_ylabel('Size (nodes)')
ax1.set_title('Heat map of egg (°C)')
fig.colorbar(im1, ax=ax1) 

# Right hand plot
x_data = np.zeros(np.size(data_000,1))
y_data = np.zeros_like(x_data)
for i in range(np.size(x_data)):
    x_data[i] = i * dx
    y_data[i] = data_000[index_reach_63, i, crosssection_point]
ax2.plot(x_data, y_data, 'b-')
ax2.set_title('Cross section')
ax2.set_xlabel('Size (mm)')
ax2.set_ylabel('Temperature (°C)')

fig.suptitle('Profile of egg at 63°C inflection point')

ax1.grid(False)
ax2.grid(False)

plt.tight_layout()
plt.style.use('seaborn-v0_8-whitegrid')
plt.show()
No description has been provided for this image

Figure 5. Heat map and cross-section of egg at 63°C inflection point.
Heat map shows 0°C boundary condition simulating egg's submersion into ice water. The profile shows the center of the yolk ≈ 30°C.

Show code
In [14]:
ani = None
count = 0
images = []

fig, (ax1, ax2) = plt.subplots(1, 2)

plt.style.use('seaborn-v0_8-whitegrid')

ax1.axhline(crosssection_point, linestyle='--', label='plot line', color='blue')
ax2.set_ylim(-5, 105)

ax1.grid(False)
ax2.grid(False)

y_label = ["boundary T", "mean white T", "mean yolk T", "min yolk T", "max yolk T", "max_white_T"]
y_colour = ["g-", "y-", "b-", "b--", "b--", "w "]

frames = np.arange(0, np.size(data_000, 0), 350)
x_data = frames * dt
y_data = np.empty((np.size(y_label), np.size(x_data)))
y_data[:, :] = np.nan

ax1.set_xlabel('Size (nodes)')
ax1.set_ylabel('Size (nodes)')
ax1.set_title('Heat map of egg (°C)')

ax2.set_title('Cross section')
ax2.set_xlabel('Time (s)')
ax2.set_ylabel('Temperature (°C)')

fig.suptitle('Profile of egg at 63°C inflection point')

for n in range(np.size(y_data, 1)):  
# for n in range(36):  # Used to generate website image

    img = ax1.imshow(data_000[frames[n]].T, cmap='hot', vmin=0, vmax=100, animated=True, origin='lower')  

    line_art = None
    # x_data[n] = frames[n] * dt
    # for i in range(np.size(y_data, 0)):
    for i in range(np.size(y_label)):
        y_data[i, n] = y_000[i, frames[n]]
        
    line_art0, = ax2.plot(x_data, y_data[0], y_colour[0], label=y_label[0], animated=True)
    line_art1, = ax2.plot(x_data, y_data[1], y_colour[1], label=y_label[1], animated=True)
    line_art2, = ax2.plot(x_data, y_data[2], y_colour[2], label=y_label[2], animated=True)
    line_art3, = ax2.plot(x_data, y_data[3], y_colour[3], label=y_label[3], animated=True)
    line_art4, = ax2.plot(x_data, y_data[4], y_colour[4], label=y_label[4], animated=True)
    fill_art   = ax2.fill_between(x_data, y_data[3], y_data[4], facecolor='grey', alpha=0.5, animated=True)

    images.append([img, line_art0, line_art1, line_art2, line_art3, line_art4, fill_art])
    count += 1
    print(f"count: {count}\r", end='')

fig.colorbar(img, ax=ax1) 
plt.tight_layout()
ani = animation.ArtistAnimation(fig, images, interval=50, blit=True, repeat_delay=2000)

plt.close() 
HTML(ani.to_jshtml())

count: 100
Out[14]:
No description has been provided for this image

Figure 6. Heat map and data from simulation.
Heat map shows 0°C boundary condition simulating egg's submersion into ice water. Green is boundary condition (temperature of water cooking the egg); Yellow is the mean temperature of the white; Blue is the mean (solid), max and min (dashed) yolk temperature. Grey area highlights temperature range of yolk. Black dash-dot lines are at 65°C and 70°C, the limits for the yolk.

Show code
In [15]:
ani = None
count = 0
images = []

fig, (ax1, ax2) = plt.subplots(1, 2)
plt.style.use('seaborn-v0_8-whitegrid')

ax1.grid(False)
ax2.grid(False)

ax1.axhline(crosssection_point, linestyle='--', label='plot line', color='blue')

ax1.set_xlabel('Size (nodes)')
ax1.set_ylabel('Size (nodes)')
ax1.set_title('Heat map of egg (°C)')

ax2.set_title('Cross section')
ax2.set_xlabel('Size (mm)')
ax2.set_ylabel('Temperature (°C)')

fig.suptitle('Profile of egg at 63°C inflection point')

for n in range(0, np.size(data_000,0), 350):  
    img = ax1.imshow(data_000[n].T, cmap='hot', vmin=0, vmax=100, animated=True, origin='lower')  

    x_data = np.zeros(np.size(data_000,1))
    y_data = np.zeros_like(x_data)
    for i in range(np.size(x_data)):
        x_data[i] = i * dx
        y_data[i] = data_000[n, i, crosssection_point]
    line_art, = ax2.plot(x_data, y_data, 'b-', animated=True)

    images.append([img, line_art])
    count += 1
    print(f"count: {count}\r", end='')

fig.colorbar(img, ax=ax1) 
plt.tight_layout()
ani = animation.ArtistAnimation(fig, images, interval=50, blit=True, repeat_delay=1000)

plt.close()
HTML(ani.to_jshtml())

count: 100
Out[15]:
No description has been provided for this image

Figure 7. Heat map and cross-section of egg at 63°C inflection point.
Heat map shows 0°C boundary condition simulating egg's submersion into ice water. The profile shows the center of the yolk.

References

[1]: https://www.sciencedirect.com/science/article/pii/S2665927124001989 - Modeling of cooking and phase change of egg white using computational fluid dynamics. R E. Sánchez-García et al. 2024

[2]: https://www.nature.com/articles/s44172-024-00334-w - Periodic cooking of eggs. E D. Lorenzo et al. 2025

[3]: https://math.arizona.edu/~gabitov/teaching/141/math_485/Final_Presentations/Cooking_Final_Presentation.pdf - Numerical Modeling of Diffusion and Phase Transitions in Heterogeneous Media. H. Jeffrey.

[4]: https://tokusen.store/en/blogs/recipes/ajitsuke-tamago-oeufs-marines-a-la-japonaise?srsltid=AfmBOop-R7FBSF8x4BBuf1P6eeOMCkhl4dBXGh3gvz_CM9kgHepp4Yd1 - Japanese Ramen Eggs (Ajitsuke Tamago). Tokusen [website]. 2025.

[4]: https://khymos.org/2009/04/09/towards-the-perfect-soft-boiled-egg/ - Towards the perfect soft boiled egg, M. Lersch. 2009

[5]: https://www.sciencedirect.com/science/article/pii/S0032579124006485 - Comparison of egg-shape equations using relative curvature measures of nonlinearity. M. Lain et al. 2024

[6]: https://newton.ex.ac.uk/teaching/CDHW/egg/ - The Science of Boiling an Egg. C. Williams.

[7]: https://www.sciencedirect.com/science/article/pii/S2665927124001989 - Modeling of cooking and phase change of egg white using computational fluid dynamics. R. Sanchez-Garcia et al. 2024

[8]: https://www.researchgate.net/figure/Thermal-properties-of-egg-components_tbl1_256804375 - Cooling of Shell Eggs with Cryogenic Carbon Dioxide: a Finite Element Analysis of Heat Transfer. C. Sabliov et al. 2002

[9]: https://www.engineeringtoolbox.com/food-thermal-conductivity-d_2177.html - The Engineering ToolBox

[10]: https://newton.ex.ac.uk/teaching/CDHW/egg/CW061201-1.pdf - Boiling an egg. C. Williams. 2006

[11]: https://douglasbaldwin.com/sous-vide.html - Dr. Douglas Baldwin, expert in sous vide cooking and applied mathematics. Perfect egg.

Appendix

Summary: Finite-difference equations

Definition of a derivative $$ \left. \frac{df(x)}{dx} \right|_{x=0} = f'(a)=\lim_{x \to 0}\frac{f(x)-f(a)}{x-a} $$

Forward difference: $$ \left. \frac{df(x)}{dx} \right|_{x_i} \approx \frac{f(x_{i+1})-f(x_i)}{x_{i+1}-x_i} $$

Backward difference: $$ \left. \frac{df(x)}{dx} \right|_{x_i} \approx \frac{f(x_i)-f(x_{i-1})}{x_i-x_{i-1}} $$

Central difference: $$ \left. \frac{df(x)}{dx} \right|_{x_i} \approx \frac{f(x_{i+1})-f(x_{i-1})}{x_{i+1}-x_{i-1}} $$

Taylor series

Assume all nodes/points are equally spaced i.e. h is constant

$$ f(x_{i+1}) = f(x_i) + hf'(x_i) + \frac{1}{2!}h^2f''(x_i) + \frac{1}{3!}h^3f'''(x_i) + \frac{1}{4!}h^4f''''(x_i) + ... $$

$$ f(x_{i+1}) = f(x_i) + hf'(x_i) + \frac{1}{2!}h^2f''(x_i) + ... + \left. \frac{1}{n!}h^n \frac{d^nf}{x^n} \right|_{x_i} + R_n \xi $$

where $\xi$ is an unknown location between $x_i \ge \xi \ge x_{i+1}$ and we cannot calculate $\xi$.

Using the first 2 terms of the Taylor's Series for $f(x_{i+1})$

$$ f(x_{i+1}) = f(x_i) + hf'(x_i) + \frac{1}{2!}h^2f''(\xi) $$

Solve for $f'(x_{i+1})$: $$ f'(x_{i+1}) = \frac{f(x_{i+1})-f(x_i)}{h} - \frac{1}{2!}hf''(\xi) = \frac{f(x_{i+1})-f(x_i)}{h} - error $$

Note: As the value of $h$ gets smaller the $error$ is smaller: $$ error = \frac{f''(\xi)}{2}h $$ This is also a first order method (order of h given to the truncation error)

General 1D diffusion problem

Diffusion type systems have a general form of: $$ \frac{\delta u}{\delta t} = \alpha \frac{\delta^2 u}{\delta x^2} + f $$

The forward Euler model1, the space domain is defined with the following nodes: $$ x_i = i\Delta x, i = 0,...,N_x $$ and the time domain is defines as: $$ t_n = n \Delta t, n = 0,...,N_t $$

So, considering the general diffusion system notes above, $u(x_i, t_n)$ can be denoted as $u^n_i$. Thus, the general diffusion system becomes: $$ \frac{\delta}{\delta t}u(x_i, t_n) = \alpha \frac{\delta^2}{\delta x^2} u(x_i, t_n) $$

The solution to a 2nd order finite difference approximation can be derived by combining 2 Taylor series, this results in a forward difference in time and a central difference in space: $$ \frac{u_i^{n+1} - u_i^n}{\Delta t} = \alpha \frac{u_{i+1}^n - 2u_i^n +u_{i-1}^n}{\delta x^2} + Err_{trunc}(dx^2) $$

Where the truncation error ($Err_{trunc}$) can be seen to be 2^nd order. Rearranging this with respect to the unknown $u_i^{n+1}$: $$ u_i^{n+1} = u_i^n + F(u_{i+1}^n - 2u_i^n + u_{i-1}^n) +\Delta t.Err_{trunc}(dx^2) $$ where: $$ F = \alpha \frac{\Delta t}{\Delta x^2} $$

F is the Fourier Number, a dimensionless value that represents the physical parameters of the diffusion problem in $\alpha$ and the discretised parameters $\Delta x$ and $\Delta t$.

[1]: Finite difference methods for diffusion processes. H. P. Langtangen, S. Linge. 2016 July 14

General 2D diffusion problem

Diffusion type systems in 2D have a general form of: $$ \frac{\delta u}{\delta t} = \alpha (\frac{\delta^2 u}{\delta x^2} + \frac{\delta^2 u}{\delta i^2}) $$ For the 2D problem we add another dimension: $$ y_j = j\Delta i, j = 0,...,N_y $$ So our 2D diffusion system, $u(x_i, y_j, t_n)$ can be denoted as $u^n_{i,j}$, thus: $$ u_i^{n+1} = u_{i,j}^n + F_x(u_{i+1,j}^n - 2u_{i,j}^n + u_{i-1,j}^n) + F_y(u_{i,j+1}^n - 2u_{i,j}^n + u_{i,j-1}^n) $$ where: $$ F_x = \alpha \frac{\Delta t}{\Delta x^2}, F_y = \alpha \frac{\Delta t}{\Delta i^2} $$ If we assume $\Delta x = \Delta i$, the 2D solution can be simplified to: $$ u_i^{n+1} = u_{i,j}^n + F(u_{i+1,j}^n + u_{i-1,j}^n + u_{i,j+1}^n + u_{i,j-1}^n - 4u_{i,j}^n) $$ where: $$ F = \alpha \frac{\Delta t}{\Delta x^2} $$

Generic finite-difference function

Show code
In [80]:
""""
Function for solving diffusion system in 2D using finite difference method.

This function has the addition of:
 - Variable alpha values for diffusion coefficient
"""
import numpy as np

def solve_diffusion_2d_vAlpha(nx, ny, dx, dy, dt, nt, alpha, u0):
    """
    Solve the 2D diffusion equation using finite difference method.

    Parameters:
    nx : int
        Number of grid points in x direction.
    ny : int
        Number of grid points in i direction.
    dx : float
        Grid spacing in x direction.
    dy : float
        Grid spacing in i direction.
    dt : float
        Time step size.
    nt : int
        Number of time steps to simulate.
    alpha : 2D array
        Diffusion coefficient.
    u0 : 2D array
        Initial condition for the diffusion equation.

    Returns:
    u : 3D array (t, x, i)
        All solutions from t = 0 to nt time steps.
    """
    # Check stability condition - using max alpha for worst case
    if dt > (dx**2 * dy**2) / (2 * alpha.max() * (dx**2 + dy**2)):
        dt_max = (dx**2 * dy**2) / (2 * alpha.max() * (dx**2 + dy**2))
        alpha_max = (dx**2 * dy**2) / (2 * dt * (dx**2 + dy**2))
        raise ValueError("Time step size dt is too large for stability. Maximum dt: ", dt_max, " or reduce alpha to: ", alpha_max)
    
    u = np.zeros((nt, nx, ny))  # Initialize the solution array
    u[0, :, :] = u0  # Set initial condition
    for n in range(1, nt):
        # Reset boundary conditions at each time step
        u[n, 0, :] = u0[0, :]       # Left boundary
        u[n, -1, :] = u0[-1, :]     # Right boundary
        u[n, :, 0] = u0[:, 0]       # Bottom boundary
        u[n, :, -1] = u0[:, -1]     # Top boundary  
        # Iterate over internal grid points
        for i in range(1, nx-1):    
            for j in range(1, ny-1):
                u[n, i, j] = (u[n-1, i, j] + 
                               alpha[i, j] * dt * ((u[n-1, i+1, j] - 2*u[n-1, i, j] + u[n-1, i-1, j]) / dx**2 +
                                              (u[n-1, i, j+1] - 2*u[n-1, i, j] + u[n-1, i, j-1]) / dy**2))
    return u

The above plot shows how the diffusion coefficient changes as temperature increases. It can be noted that:

  • The shell had a constant thermal diffusion coefficient
  • The white's thermal diffusion coefficient decreases at the 63℃ transition temperature
    • Note: This is assumed to be the transition due to the denaturing of the ovotransferrin protein as the later denaturing is not until 80℃.
    • Therefore further change in diffusion coefficient > 80℃ are not covered in the data for this model.
  • The yolk's thermal diffusion coefficient increases at the 70℃ transition temperature.
Show code
In [83]:
""""
Function for solving diffusion system in 2D using finite difference method.

This function has the addition of:
 - Variable alpha values for diffusion coefficient
"""
import numpy as np

def solve_diffusion_2d_vAlpha(nx, ny, dx, dy, dt, nt, alpha, u0):
    """
    Solve the 2D diffusion equation using finite difference method.

    Parameters:
    nx : int
        Number of grid points in x direction.
    ny : int
        Number of grid points in i direction.
    dx : float
        Grid spacing in x direction.
    dy : float
        Grid spacing in i direction.
    dt : float
        Time step size.
    nt : int
        Number of time steps to simulate.
    alpha : 2D array
        Diffusion coefficient.
    u0 : 2D array
        Initial condition for the diffusion equation.

    Returns:
    u : 3D array (t, x, i)
        All solutions from t = 0 to nt time steps.
    """
    # Check stability condition - using max alpha for worst case
    if dt > (dx**2 * dy**2) / (2 * alpha.max() * (dx**2 + dy**2)):
        dt_max = (dx**2 * dy**2) / (2 * alpha.max() * (dx**2 + dy**2))
        alpha_max = (dx**2 * dy**2) / (2 * dt * (dx**2 + dy**2))
        raise ValueError("Time step size dt is too large for stability. Maximum dt: ", dt_max, " or reduce alpha to: ", alpha_max)
    
    u = np.zeros((nt, nx, ny))  # Initialize the solution array
    u[0, :, :] = u0  # Set initial condition
    for n in range(1, nt):
        # Reset boundary conditions at each time step
        u[n, 0, :] = u0[0, :]       # Left boundary
        u[n, -1, :] = u0[-1, :]     # Right boundary
        u[n, :, 0] = u0[:, 0]       # Bottom boundary
        u[n, :, -1] = u0[:, -1]     # Top boundary  
        # Iterate over internal grid points
        for i in range(1, nx-1):    
            for j in range(1, ny-1):
                u[n, i, j] = (u[n-1, i, j] + 
                               alpha[i, j] * dt * ((u[n-1, i+1, j] - 2*u[n-1, i, j] + u[n-1, i-1, j]) / dx**2 +
                                              (u[n-1, i, j+1] - 2*u[n-1, i, j] + u[n-1, i, j-1]) / dy**2))
    return u